Prévision de la performance énergétique des bâtiments


La performance énergétique est un enjeu majeur de la réduction de la consommation d’énergie et des émissions de gaz à effet de serre. La menace du changement climatique pèse sur l’humanité et la transition énergétique vise notamment à améliorer cette performance énergétique des bâtiments. Depuis plus d’une décennie, le diagnostic de performance énergétique vient faciliter ce travail en apportant des données clés.

L'étude menée se distingue en 2 phases :

  • Exploration statistique pour comprendre les différentes variables explicatives du jeu de données et observer des éventuelles corrélations.
  • Modélisation des performances énergétiques de nouveaux batiments à l'aide d'algorithmes prédictifs (Random Forest, SVM, etc). Il est important de noter que nous allons traiter le problème selon 2 aspects : Classification et Régression.


Importation des packages

In [113]:
library(corrplot)
library(ggplot2)
library(lattice)
library(rpart)
library(partykit)
library(randomForest)
library(ROCR)
library(missForest)
library(reshape2)
library(FactoMineR)
library(factoextra)
#library(rgl)
library(plotly)
library(dplyr)
library(gridExtra)
library(bestglm)
library(MASS)
library(VGAM)
library(e1071)
library(glmnet)
library(broom)
library(gbm)
library(bestglm)
library(adabag)
library(caret)

Importation des données

In [2]:
data.complet <- read.csv('DataEnergy.csv')
x <- data.complet
head(x)
Relative.compactnessSurface.areaWall.areaRoof.areaOverall.heightorientationGlazing.areaGlazing.area.distrEnergyEnergy.efficiency
0.9829276 530.4900 306.4846 112.0027 7 North 1.609490e-020 34.26394 B
0.9835473 519.8724 299.7763 110.0480 7 East -9.386813e-030 34.58975 B
0.9794535 516.1912 303.3744 106.4084 7 South -6.974937e-040 38.77805 C
0.9777325 518.9241 292.8122 113.0559 7 West 9.554434e-060 37.94781 C
0.9030294 552.9689 316.2361 118.3664 7 North -6.592326e-030 47.67586 D
0.8909102 558.6037 314.9162 121.8437 7 East -1.612145e-020 41.90847 C

Descriptif des variables

  • Relative compactness : compacité relative (propriété liée à la forme du batiment et sa surface).
  • Surface area : surface totale du batiment.
  • Wall area : surface des murs.
  • Roof area : surface du toit.
  • Overal height : hauteur du batiment (variable qualitative ordonnée: 3,5m et 7m).
  • Orientation : orientation du batiment (North, East, South, West).
  • Glazing area : surface totale des vitrages du batiment.
  • Glazing area distr : orientation des vitrages du batiment (55% East, 55% South, etc).
  • Energy : quantifie les performances énergétique du batiment.
  • Energy efficiency : lettre de classification des performances énergétique du batiment.

Préparation des variables

In [3]:
str(x)
'data.frame':	768 obs. of  10 variables:
 $ Relative.compactness: num  0.983 0.984 0.979 0.978 0.903 ...
 $ Surface.area        : num  530 520 516 519 553 ...
 $ Wall.area           : num  306 300 303 293 316 ...
 $ Roof.area           : num  112 110 106 113 118 ...
 $ Overall.height      : num  7 7 7 7 7 7 7 7 7 7 ...
 $ orientation         : Factor w/ 4 levels "East","North",..: 2 1 3 4 2 1 3 4 2 1 ...
 $ Glazing.area        : num  1.61e-02 -9.39e-03 -6.97e-04 9.55e-06 -6.59e-03 ...
 $ Glazing.area.distr  : int  0 0 0 0 0 0 0 0 0 0 ...
 $ Energy              : num  34.3 34.6 38.8 37.9 47.7 ...
 $ Energy.efficiency   : Factor w/ 7 levels "A","B","C","D",..: 2 2 3 3 4 3 3 3 4 4 ...
In [4]:
# Transformation des variables qualitatives en variables catégorielles


x$Glazing.area.distr <- as.factor(x$Glazing.area.distr)
levels(x$Glazing.area.distr) <- c("No Glazing", "55% North", "55% East", 
                                  "55% South", "55% West", "Uniform")


x$Energy.efficiency <- factor(x$Energy.efficiency, ordered=TRUE)
In [5]:
summary(x)
 Relative.compactness  Surface.area     Wall.area       Roof.area    
 Min.   :0.6125       Min.   :501.4   Min.   :234.3   Min.   :105.3  
 1st Qu.:0.6779       1st Qu.:598.7   1st Qu.:291.8   1st Qu.:137.4  
 Median :0.7517       Median :673.1   Median :315.8   Median :183.3  
 Mean   :0.7645       Mean   :671.3   Mean   :318.3   Mean   :176.5  
 3rd Qu.:0.8350       3rd Qu.:744.6   3rd Qu.:343.0   3rd Qu.:220.5  
 Max.   :0.9912       Max.   :826.0   Max.   :425.8   Max.   :225.8  
                                                                     
 Overall.height orientation  Glazing.area       Glazing.area.distr
 Min.   :3.50   East :192   Min.   :-0.02401   No Glazing: 48     
 1st Qu.:3.50   North:192   1st Qu.: 0.10306   55% North :144     
 Median :5.25   South:192   Median : 0.24752   55% East  :144     
 Mean   :5.25   West :192   Mean   : 0.23435   55% South :144     
 3rd Qu.:7.00               3rd Qu.: 0.39119   55% West  :144     
 Max.   :7.00               Max.   : 0.42702   Uniform   :144     
                                                                  
     Energy      Energy.efficiency
 Min.   :10.21   A:208            
 1st Qu.:29.36   B:109            
 Median :41.76   C: 80            
 Mean   :46.92   D: 79            
 3rd Qu.:64.33   E:109            
 Max.   :94.84   F:102            
                 G: 81            

On remarque des valeurs négatives pour la variable Glazing Area, qui représente pourtant la surface totale des vitrages. En fait, ces valeurs négatives sont dues à un bruit ajouté. Regardons comment sont réparties ces valeurs négatives afin de leur appliquer une correction.


In [6]:
summary(x[which(x$Glazing.area<0),]$Glazing.area.distr)
No Glazing
24
55% North
0
55% East
0
55% South
0
55% West
0
Uniform
0

D'après le resultat ci-dessus, toutes les valeurs de Glazing.area inférieures à 0 appartiennent à la catégorie No glazing. On va donc forcer leur valeur à 0.


In [7]:
x$Glazing.area[x$Glazing.area < 0] <- 0


On remarque aussi que la variable Overall.height prend seulement deux niveaux ordonnés : 3,5m et 7m. On va donc considérer cette variable comme catégorielle ordonnée.

In [8]:
x$Overall.height <- as.ordered(x$Overall.height)


A. Analyse des données

1. Statistiques descriptives univariées


In [9]:
options(repr.plot.width = 4, repr.plot.height = 4) 
histogram(x$Energy.efficiency, freq=FALSE, xlab="Energy Efficiency", 
          main="Proportion du nombre de batiments \npar classe énergétique")

La répartition des données par classe énergie est plutot uniforme. On observe cependant que la classe A est majoritaire dans ce jeu de données.


In [10]:
options(repr.plot.width = 6, repr.plot.height = 6) 
boxplot(x[,-c(5,6,8,10,11)], col='lightblue', xaxt = "n", 
        main="Box plot des variables quantitatives") 
grid()
text(seq_along(x[,-c(5,6,8,10,11)]), par("usr")[3] - 0.5, labels = names(x[,-c(5,6,8,10,11)]), srt = 45, adj = 1, xpd = TRUE);

Les variables ne sont pas toutes de même ordre de grandeur mais les distributions semblent homogènes. Etant donné les ordres de grandeurs bien différents, une attention particulière sera accordée lors de la réalisation de l'ACP.

Il faut maintenant étudier plus précisemment la distribution de ces variables.


In [11]:
par(mfrow=c(2,3))
options(repr.plot.width = 9, repr.plot.height = 6) 

hist(x$Glazing.area, xlab='Glazing Area', main='Glazing Area Distribution',col='lightblue')
hist(x$Relative.compactness,xlab= 'Relative Compactness' ,main = 'Relative Compactness Distribution',col='lightblue')
hist(x$Roof.area,xlab= 'Roof Area' ,main = 'Roof Area Distribution',col='lightblue')
hist(x$Surface.area,xlab= 'Surface Area' ,main = 'Surface Distribution',col='lightblue')
hist(x$Wall.area,xlab= 'Wall Area' ,main = 'Wall Area Distribution',col='lightblue')
hist(x$Energy,xlab= 'Energy' ,main = 'Energy Distribution',col='lightblue')

Il n'y a pas de transformation de variables nécessaires. On remarque tout de même que la variable Roof Area est inéquitablement répartie.


In [12]:
options(repr.plot.width = 4, repr.plot.height = 4) 
pie(table(x$orientation), main = "Proportion du nombre de batiment \n selon l'orientation")

La répartition des orientations des bâtiments est uniforme pour le jeu de données.


In [13]:
options(repr.plot.width = 6, repr.plot.height = 3) 
histogram(x$Glazing.area.distr,xlab= 'Glazing area distr' ,
          main = 'Proportion du nombre de batiments \n par orientation des vitrages',col='gold')

La répartition des orientations des vitrages est équitable entre les différentes classes (55% Nord, 55% Est, 55% Sud, 55% Ouest et sans vitrage). Cependant, les batiments sans vitres sont minoritaires.


2. Statistiques descriptives multivariées

In [14]:
options(repr.plot.width = 6, repr.plot.height = 5) 
M <- cor(x[,-c(5,6,8,10,11)])
corrplot(M)
  • Relative compactness est très fortement correlé négativement avec Surface area (-0.98) et Roof Area (-0.87). Il doit exister un lien numérique entre ces variables que nous étudierons par la suite.

  • Roof Area est aussi fortement corrélé avec Energy.

Il est difficile de tirer des conclusions avec le graphique des corrélations. Regardons la matrice de scatterplots.


In [15]:
options(repr.plot.width = 6, repr.plot.height = 6)
pairs(x[,-c(5,6,8,10,11)])


On remarque une relation linéaire entre Relative.compactness et Surface.area.

On voit que Glazing.area, qui est une variable quantitative, est finalement répartie selon 4 intervalles.


In [16]:
defaultW <- getOption("warn")
options(warn = -1)

fig <- plot_ly(x, x = ~Roof.area, y = ~Wall.area, z = ~Surface.area, color = ~Energy.efficiency)
fig <- fig %>% add_markers()
fig <- fig %>% layout(scene = list(xaxis = list(title = 'Roof Area'),
                     yaxis = list(title = 'Wall Area'),
                     zaxis = list(title = 'Surface Area')))

fig
options(warn = defaultW)
In [17]:
## On vérifie la relation Surface.area = Wall.area + 2*Roof.area

a=x$Surface.area-(x$Wall.area+2*x$Roof.area)
cat("Nombre de valeurs pour lesquelles la relation 'Surface.area = Wall.area + 2*Roof.area' n'est pas vérifié : ", length(which(round(a,3)!=0)))
Nombre de valeurs pour lesquelles la relation 'Surface.area = Wall.area + 2*Roof.area' n'est pas vérifié :  0


La surface est calculée en sommant la surface des murs et en ajoutant deux fois la surface du toit. La relation étant vérifiée pour tous les points de notre jeu de donnée, nous allons pouvoir supprimer une ou plusieurs variables étant donné que les informations sont redondantes.


In [18]:
options(repr.plot.width = 18, repr.plot.height = 7) 

par(mfrow=c(1,3))
boxplot(Wall.area ~Energy.efficiency, data=x, col='lightblue', main="Wall Area vs Energy Efficiency")
boxplot(Roof.area ~Energy.efficiency, data=x, col='lightgreen', main="Roof Area vs Energy Efficiency")
boxplot(Surface.area ~Energy.efficiency, data=x, col='lightyellow', main="Surface Area vs Energy Efficiency")

On remarque sur ces graphiques que les classes A, B et C ont une surface de toît élevée. De la même façon, ces classes ont une surface totale élevée.

La tendance est inversée pour les classes D, E, F et G.

Pour continuer l'analyse, on va donc supprimer la variable Roof Area, d'une part car son information est contenue dans la variable Surface Area et d'autre part car sa distribution est inéquitablement répartie.

On retire donc la variable Roof Area.

In [19]:
x <- dplyr::select(x,-Roof.area)
In [20]:
options(repr.plot.width = 10, repr.plot.height = 5) 

par(mfrow=c(1,2))
boxplot(x$Overall.height~ x$Energy.efficiency, main ="Hauteur \n selon la classe Energetique", xlab="Classe Energetique", ylab="Hauteur",col='red')
boxplot(x$Relative.compactness~ x$Energy.efficiency, main ="Relative compactness \n selon la classe energetique", xlab="Classe Energetique", ylab="Hauteur",col='darkred')

On remarque sur le graphique que les classes A, B et C ont majoritairement une hauteur de 3m50 tandis que les classes D, E, F et G ont principalement une hauteur plus élevée (7m).


De la même façon, on remarque que la valeur de Relative Compactness est plus faible pour les bonnes classes énergetiques par rapport aux moins bonnes.

In [21]:
options(repr.plot.width = 12, repr.plot.height = 5) 
h1 = histogram(~Energy.efficiency | orientation , data=x, xlab= "Energy Efficiency",main="Orientation vs Class Energy",)
h2 = histogram(~Energy.efficiency | Glazing.area.distr, data=x,main="Glazing area distr vs Class Energy", xlab= "Energy Efficiency")
grid.arrange(h1, h2, nrow=1, ncol=2)
#histogram(~orientation | Energy.efficiency , data=x, xlab= "Energy Efficiency",main="Orientation vs Class Energy",)
#mosaicplot(table2, main="Glazing area distr vs Class Energy")

Les classes d'énergie sont uniformement réparties selon l'orientation du bâtiment et selon l'orientation des vitrages. L'orientation des bâtiments n'influe donc probablement pas l'appartenance à une certaine classe d'énergie.

En revanche, lorsqu'il y a pas de vitrage, la classe A est clairement majoritaire. On peut faire le lien avec la réalité où on suppose qu'un batiment sans vitrage est mieux isolé.


3. Analyse en Composantes Principales

On enlève la variable Energy car c'est celle que l'on veut prédire. On garde seulement les variables quantitatives.

Remarque : Les variables n'étant pas dans les mêmes unités, nous faisons une ACP centrée réduite.

In [22]:
X_acp <- dplyr::select(x,c(Relative.compactness, Surface.area, Wall.area, Glazing.area,Energy.efficiency))
res.pca <- PCA(X_acp, scale.unit=TRUE, quali.sup=5, graph = FALSE)
In [23]:
#fviz_eig(res.pca, addlabels = TRUE, ylim = c(0, 50), new.plot=FALSE, graph.type="classic")
lbls = paste("Dimension",1:4, ":", round(res.pca$eig[,2],2), "%" )
coul = c('palegreen','salmon','lightyellow','lightblue','black')
options(repr.plot.width = 8, repr.plot.height = 9) 
par(mfrow=c(2,1))
barplot(res.pca$eig[,2],names.arg=paste("Dim.",1:nrow(res.pca$eig)), xlab="Dimensions", ylab = "% of Variance", col = coul,main="Pourcentage de variance expliquée par dimension")
legend(x='topright',legend=lbls, fill= coul, bty='n', cex=0.8, text.font=0.5)
boxplot(res.pca$ind$coord, outlier=TRUE, col=coul, main = "Distribution des premières composantes")

Au vu du pourcentage de variance expliqué et des ditributions des trois premières dimensions, nous sommes tentés de les garder.


In [24]:
options(repr.plot.width = 15, repr.plot.height = 6) 
plot1 <- plot(res.pca,choix="ind", label="none", new.plot=FALSE, graph.type="ggplot",habillage=5)
plot2 <- plot(res.pca,choix="var", new.plot=FALSE, graph.type="ggplot")
grid.arrange(plot1, plot2, nrow=1, ncol=2)

On observe des clusters dans le graphe des individus qui ne reflètent pas exactement les classes. Cependant on peut séparer les "bonnes classes" (A, B et C) des moins bonnes classes (D,E,F et G) par la première composante principale. En effet, les "bonnes classes" semblent se distinguer des "mauvaises classes", par une valeur positive sur l'axe 1. Ce dernier étant porté par Surface Area positivement et par Relative Compactness négativement.

La dimension 2 est principalement portée par Glazing Area. On ne distingue pas les différentes classes grâce à cette composante (et donc à cette variable).


In [25]:
options(repr.plot.width = 15, repr.plot.height = 6) 
plot3 <- plot(res.pca,axes=c(2,3),choix="ind", label="none", new.plot=FALSE, graph.type="ggplot",habillage=5)
plot4 <- plot(res.pca,axes=c(2,3),choix="var", new.plot=FALSE, graph.type="ggplot")
grid.arrange(plot3, plot4, nrow=1, ncol=2)

On remarque encore ici des clusters selon l'axe 2; ces derniers sont sans doute dus à la variable Glazing Area qui est disparate comme nous l'avons souligné précedemment.

On distingue tout de même que les bonnes classes ont tendance à avoir une valeur négative sur l'axe 3, porté par Wall Area. La variable Wall Area permet sans doute elle aussi de distinguer certaines classes, plus Wall Area est faible, plus la classe Energétique est bonne.

Les conclusions de l'ACP correspondent aux caractéristiques montrées précedemment sur l'impact de certaines variables sur la consomation Energetique.


4. Clustering

In [26]:
x_num=dplyr::select(x, Relative.compactness, Wall.area,Surface.area, Glazing.area)
x_num_cr=scale(x_num)
In [27]:
stacked_cluster = hclust(dist(x_num_cr), method = "ward.D")
plot(stacked_cluster)
In [28]:
# Elbow method
fviz_nbclust(x_num_cr, kmeans, method = "wss") +
    geom_vline(xintercept = 2, linetype = 2)+
  labs(subtitle = "Elbow method")

Ces deux méthodes nous suggèrent 2 clusters, or on aimerait qu'il y en ai 7 (correspondant aux 7 classes énergétiques)

In [29]:
set.seed(10)
km.out2 = kmeans(x_num_cr,centers=2)
km.out7 = kmeans(x_num_cr,centers=7)
plot.kmean2<-fviz_cluster(km.out2, x_num_cr, ellipse.type = "norm", main='Kmeans with 2 clusters')
plot.kmean7<-fviz_cluster(km.out7, x_num_cr, ellipse.type = "norm", main='Kmeans with 7 clusters')
grid.arrange(plot.kmean2, plot.kmean7, nrow=1, ncol=2)
In [30]:
cat('Table de contingence avec 2 clusters :')

print(table(km.out2$cluster,x$Energy.efficiency))

cat(' \n\n ')
cat('Table de contingence avec 7 clusters :')
print(table(km.out7$cluster,x$Energy.efficiency))
Table de contingence avec 2 clusters :   
      A   B   C   D   E   F   G
  1   0   4  11  74 101  77  53
  2 208 105  69   5   8  25  28
 

 Table de contingence avec 7 clusters :   
      A   B   C   D   E   F   G
  1  69  45  38   1   0   0   0
  2   0   0   0   5  49  55  29
  3   0   0   0   3   8  25  28
  4 104  16   0   0   0   0   0
  5   0   2   4  31  11  16  24
  6   0   2   7  38  41   6   0
  7  35  44  31   1   0   0   0

Kmeans ne semble pas adapté à notre problème. Cependant, cet algorithme suggère de regrouper les classes A,B et C ensemble ainsi que D,E,F,G (bonne consommation et mauvaise consommation). Nous nous intéresserons sur cela lors d'une partie Bonus.


B. Séparation des données :

On sépare les données en un ensemble d'apprentissage (train) et un ensemble test. Le modèle est entrainé sur l'échantillon train et il est évalué sur l'échantillon test.

Cette étape est nécessaire pour permettre d'évaluer le modèle déterminé avec les échantillons d'entrainement. En effet, on peut par exemple détecter un overfitting grâce à l'échantillon test (un très bon score avec l'échantillon d'entrainement ne signifie pas toujours que c'est un très bon modèle).

In [31]:
## 75% of the sample size
x_fin<-x
train_ratio=0.7
smp_size <- floor(train_ratio * nrow(x))

## set the seed to make your partition reproducible
set.seed(2512)
train_ind <- sample(seq_len(nrow(x)), size = smp_size)

train <- x[train_ind, ]
test <- x[-train_ind, ]
In [32]:
cat("Nombre d'individus dans l'échantillon train : ", nrow(train), '\n')
cat("Nombre d'individus dans l'échantillon test : ", nrow(test))
Nombre d'individus dans l'échantillon train :  537 
Nombre d'individus dans l'échantillon test :  231
In [33]:
options(repr.plot.width = 5, repr.plot.height = 4) 
histogram(train$Energy.efficiency, freq=FALSE, xlab="Energy Efficiency", 
          main="Proportion du nombre de batiments \npar classe énergétique dans le train set")

On remarque une plus forte proportion d'individus de la classe A. Nos modèles auront donc tendance à être meilleurs pour la classe A.

Nous avons choisi des métriques de performance qui tiennent compte de ce fait.


C. Modélisation

Définition des métriques

Pour évaluer les performances d'un modèle de classification, il faut définir des mesures traduisant son efficacité:

Quatres scores ont été utilisées pour comparer les différents modèles :

  • Score accuracy : Quotient entre les bien prédits et le nombre total de prédictions (fonction perf_model). Le défaut de cette première métrique est qu'elle ne prend pas en compte l'écart de classe quand on se trompe. En effet, imaginons que la vraie valeur soit la classe A ; l'erreur est bien plus grave si on prédit une consommation de type G que si l'on prédit une consommation de type B.

  • Root Mean Squared Error (RMSE) : Uniquement utilisé dans le cadre de la régression

  • Mean Absolute Error (MAE) : Cette métrique est obtenue en utilisant la formule trouvée dans le papier Measuring the performance of ordinal classification (Aime S. Cardoso & Ricardo Sousa). Bien qu'associé en temps normal à des regréssions, cette mesure est ici adaptée pour prendre en compte les écart de classe lors d'une classification. (fonction MAE()). Plus le score est faible, meilleur est la classification.

  • Macro F1 calcule la métrique F1 pour chaque classe individuellement (comme pour une classification binaire où la valeur TRUE correspondrait à la classe que l'on souhaite prédire la valeur FALSE regrouperait toutes les autres classes.). On fait ensuite la moyenne de toutes ces mesures F1. Il est important de noter que cette métrique ne tient pas en compte des classes non uniformément réparties.

NB : Nous gardons aussi toutes les valeurs individuelles de mesures F1 et nous les afficherons à la fin pour comparer nos modèles. Certaines méthodes (linéaire ou non-linéaires) seront en effet sans doutes meilleures pour prédire certaines classes que d'autres.

Toutes les métriques pour la multiclassification peuvent être trouvées ici :

http://www.inescporto.pt/~jsc/publications/journals/2011JaimeIJPRAI.pdf

In [34]:
## Définition des scores

perf_model <- function(table_contingence){
    return( sum(diag(table_contingence))/sum(table_contingence))
}

MAE <- function(table){
    lignes = nrow(table)
    cols = ncol(table)
    N = sum(table)
    MSE=0
    for (i in seq(1,lignes)){
        
        for (j in seq(1,cols)){
            
            MSE = MSE + (table[i,j]*abs(i-j))
            
        }
    }
    
 return(MSE/N)   
}

macro_F1 <- function(table){
    lignes = nrow(table)
    cols = ncol(table)
    diag = diag(table)
    rowsums = apply(table, 1, sum)   
    colsums = apply(table, 2, sum) 
    precision = diag / colsums 
    recall = diag / rowsums 
    f1 = 2 * precision * recall / (precision + recall)   
    
    return(f1)
}
In [35]:
#Definition des tableaux de score :
#Ces fonctions permettent de regrouper dans un tableau les différents scores pour chaque modèle
Tab_score_class = as.data.frame(setNames(replicate(4,numeric(0), simplify = F),
                                         c("Modele","Accuracy","MAE", "Macro F1") ))
Tab_score_reg = as.data.frame(setNames(replicate(5,numeric(0), simplify = F),
                                       c("Modele","Accuracy","MAE", "Macro F1","RMSE") ))

Tab_F1_class = as.data.frame(setNames(replicate(8,numeric(0), simplify = F),
                                         c("Modele","A","B","C","D","E","F","G") ))

Tab_F1_reg = as.data.frame(setNames(replicate(8,numeric(0), simplify = F),
                                         c("Modele","A","B","C","D","E","F","G") ))

#pred est la liste des prédictions sur le jeu de test
#quali est le type de la prediction

Compute_Error <- function(pred,quali=TRUE, name_model=""){ 
    
    if (quali){
        table_result = table(pred.reg = pred , observations = test$Energy.efficiency)
        cat("Table de contingence de", name_model, ": \n")
        print(table_result)
        score1 = perf_model(table_result)
        score2 = MAE(table_result)
        score3 = mean( macro_F1(table_result) )
        #on cree une df avec les mêmes noms de colonnes
        to_add <- data.frame(Modele=name_model, Accuracy=score1,MAE=score2, MacroF1=score3)
        Tab_score_class = rbind(Tab_score_class,to_add) # on concatène ensuite
        return(Tab_score_class)
    }
    else{
        pred_class <- classify(pred)   
        table_result = table(pred.reg = pred_class , observations = test$Energy.efficiency)
        cat("Table de contingence de", name_model, ": \n")
        print(table_result)
        score1 = perf_model(table_result)
        score2 = MAE(table_result)
        score3 = mean( macro_F1(table_result) )
        rmse = sqrt( sum( (pred - test$Energy)^2 ) /nrow(test) )
        #on cree une df avec les mêmes noms de colonnes
        to_add <- data.frame(Modele=name_model, Accuracy=score1,MAE=score2, MacroF1=score3, RMSE=rmse) 
        Tab_score_reg = rbind(Tab_score_reg,to_add) # on concatène ensuite
        return(Tab_score_reg)        
        
    }
}


Compute_F1 <- function(predict,quali=TRUE, name_model=""){
    
    if (quali){
        table_result = table(pred.reg = predict , observations = test$Energy.efficiency)
        temp = macro_F1(table_result)
        add_F1= data.frame(Modele = name_model,"A"= temp[[1]],"B"= temp[[2]], "C"= temp[[3]], "D"=temp[[4]],
                           "E"=temp[[5]],"F"=temp[[6]],"G"=temp[[7]])
        
        
        Tab_F1_class = rbind(Tab_F1_class,add_F1)
        return(Tab_F1_class)

    }    
    else{
        pred_class <- classify(predict)   
        table_result = table(pred.reg = pred_class , observations = test$Energy.efficiency)
        temp = macro_F1(table_result)
        add_F1= data.frame(Modele=name_model,"A"=temp[1],"B"=temp[2], "C"=temp[3], "D"=temp[4],
                           "E"=temp[5],"F"=temp[6],"G"=temp[7])
        
        Tab_F1_reg = rbind(Tab_F1_reg,add_F1)
        return(Tab_F1_reg)
    }   
}

1. Classification

Dans cette partie, nous nous concentrons uniquement sur le problème de classification, qui consiste à prédire la classe énergétique d'un batiment. La variable quantitative Energy est donc supprimée de cette partie.

In [36]:
train_class <- dplyr::select(train,-Energy)
test_class <- dplyr::select(test, -Energy)
In [37]:
x_train_class<-dplyr::select(train_class, -Energy.efficiency)
y_train_class<- train_class$Energy.efficiency

x_test_class<-dplyr::select(test_class,-Energy.efficiency)
y_test_class<-test_class$Energy.efficiency

a. Régression logistique polythomiale :

Nous allons réaliser une régression logistique polynomial ordonnée aussi appelée "polythomial regression".

Il y a deux moyens de le réaliser : additive logits and adjacents logits.

Mise en place du modèle:

In [38]:
#Additive Simple: 
vglm.c <- vglm(Energy.efficiency ~., data=train_class, family=cumulative(parallel=T, reverse=F))

Résultats :

In [39]:
p <- predict(vglm.c, newdata = test_class, type = "response") #matrice de proba d'appartenance à chaque classe
vglm.c.pred <- apply(p,1, which.max) 


Tab_score_class <-Compute_Error(vglm.c.pred,quali=TRUE,name_model = "Naive Logit");Tab_score_class
Tab_F1_class <- Compute_F1(vglm.c.pred,quali=TRUE,name_model = "Naive Logit")
Table de contingence de Naive Logit : 
        observations
pred.reg  A  B  C  D  E  F  G
       1 54 18  5  0  0  0  0
       2  2 11  8  0  0  0  0
       3  5  8  6  1  0  0  0
       4  0  2  3  8  6  0  0
       5  0  0  0  7 13 13  3
       6  0  0  0  2  9 13 10
       7  0  0  0  2  0  7 15
ModeleAccuracyMAEMacroF1
Naive Logit0.5194805 0.5714286 0.459497

Commentaire :

Le score de classification est relativement faible, mais il n'est pas mauvais. On remarque par ailleurs que les erreurs de classification sont majoritairement dans les classes voisines. Ce résultat n'est pas choquant et c'est un bon début. Dans la partie suivante, nous allons essayer de pénaliser cette régression.


b. Regression logistique pénalisée :

Mise en place du modèle :

In [40]:
x.mat.train <- model.matrix(Energy.efficiency ~ . -1 , data = train_class) #permet de gérer les variables catégorielles
lasso.c.cv <- cv.glmnet(y = y_train_class, x = x.mat.train,family="multinomial")

plot(lasso.c.cv)
paste("Estimation CV de lambda : ", round(lasso.c.cv$lambda.min,3))
'Estimation CV de lambda : 0.001'

Résultats :

In [41]:
x.mat.test <- model.matrix(Energy.efficiency ~ . -1, data = test_class)
lasso.c.cv.pred <- predict(lasso.c.cv,s="lambda.min", newx = x.mat.test,type="class") #valeur predite sur le jeu de train
Tab_score_class <-Compute_Error(lasso.c.cv.pred,quali=TRUE,name_model = "Penalised Logit");Tab_score_class
Tab_F1_class <- Compute_F1(lasso.c.cv.pred,quali=TRUE,name_model = "Penalised Logit")
Table de contingence de Penalised Logit : 
        observations
pred.reg  A  B  C  D  E  F  G
       A 52 21  5  0  0  0  0
       B  7  7  6  0  0  0  0
       C  2  8  8  0  0  0  0
       D  0  3  3 10  7  4  1
       E  0  0  0  6 10  6  6
       F  0  0  0  2  7 13  5
       G  0  0  0  2  4 10 16
ModeleAccuracyMAEMacroF1
Naive Logit 0.5194805 0.5714286 0.4594970
Penalised Logit0.5021645 0.6363636 0.4465665

Commentaires :

La pénalisation LASSO permet de diminuer le nombre de variables utilisées tout en dégradant très peu l'information. Dans ce modèle une contrainte est ajoutée sur les paramètres. La pénalisation LASSO a donc une moins bonne accuracy que le modèle linéaire car le modèle linéaire minimise uniquement les moindres carrés (et ne pénalise pas les variables).

Nous avons donc un moins bon modèle mais ce modèle est plus parcimonieux : il contient moins de variables.

Essayons maintenant des méthodes non-linéaires.


c. Arbre de décision

Mise en place du modèle optimisé

In [42]:
options(repr.plot.width = 15, repr.plot.height = 10)
tree.c=rpart(Energy.efficiency~.,data=train_class,control=rpart.control(cp=0.001))
xmat=xpred.rpart(tree.c, xval = 10)
#levels(datappr[,"Survived"]) <- c(1,2)
xerr <- (xmat-as.numeric(y_train_class))^2
CVerr=apply(xerr,2,sum)
cp_opti = as.numeric(attributes(which.min(CVerr))$names); 
cat('cp opti :', cp_opti)

tree.c.opti=rpart(Energy.efficiency~.,data=train_class,control=rpart.control(cp=cp_opti))
options(repr.plot.width = 15, repr.plot.height = 10)
plot(tree.c.opti)
text(tree.c.opti)
cp opti : 0.001601282

Commentaire : On remarque que la variable Overall.height permet de classifier les extrêmes (classe A,B vs D,E,F,G), la délimitation dela classe C n'étant pas aussi nette.

Résultats :

In [43]:
tree.c.opti.reg <- predict(tree.c.opti, newdata = test_class, type = "class")
Tab_score_class <-Compute_Error(tree.c.opti.reg,quali=TRUE,name_model = "Tree opti");Tab_score_class
Tab_F1_class <- Compute_F1(tree.c.opti.reg,quali=TRUE,name_model = "Tree Opti")
Table de contingence de Tree opti : 
        observations
pred.reg  A  B  C  D  E  F  G
       A 48 10  1  0  0  0  0
       B  8  7  4  0  0  0  0
       C  5 20 15  0  0  0  0
       D  0  2  2 15  2  1  0
       E  0  0  0  4 12  5  0
       F  0  0  0  1 14 18  7
       G  0  0  0  0  0  9 21
ModeleAccuracyMAEMacroF1
Naive Logit 0.5194805 0.5714286 0.4594970
Penalised Logit0.5021645 0.6363636 0.4465665
Tree opti 0.5887446 0.4545455 0.5638029

Commentaire :

L'apport de cette méthode est significative car notre prédiction s'améliore pour chacune des métriques étudiées. L'utilisation de méthodes non-linéaires semble pour l'instant être une bonne piste. Nous allons donc par la suite en étudier d'autres. On se dirige donc vers un RandomForest qui devrait prodiguer d'encore meilleurs résultats.


d. Random Forest

Mise en place du modèle :

In [44]:
rf.c<-randomForest::randomForest(x = x_train_class, y = y_train_class,
                         xtest = x_test_class, ytest = y_test_class,
                         importance = TRUE)

Résultats :

In [45]:
rf.c.pred <- rf.c$test$predicted
Tab_score_class <-Compute_Error(rf.c.pred,quali=TRUE,name_model = "Random Forest");Tab_score_class
Tab_F1_class <- Compute_F1(rf.c.pred,quali=TRUE,name_model = "RF");
Table de contingence de Random Forest : 
        observations
pred.reg  A  B  C  D  E  F  G
       A 55 20  3  0  0  0  0
       B  6 11  6  0  0  0  0
       C  0  6 11  0  0  0  0
       D  0  2  2 14  2  0  0
       E  0  0  0  6 13  7  0
       F  0  0  0  0 13 16  4
       G  0  0  0  0  0 10 24
ModeleAccuracyMAEMacroF1
Naive Logit 0.5194805 0.5714286 0.4594970
Penalised Logit0.5021645 0.6363636 0.4465665
Tree opti 0.5887446 0.4545455 0.5638029
Random Forest 0.6233766 0.3982684 0.5929760

Commentaire : Nous remarquons dans un premier temps que le taux de classification est meilleur que celui obtenu avec l'arbre de décision.

Pour améliorer le taux de classification de Random Forest, nous pouvons essayer de modifier le paramètre mtry correspondant au nombre de variables testées à chaque division.

Utilisons la librairie caret.

Mise en place du modèle optimisé :

In [46]:
cvControl <- trainControl(method = "cv", number = 10)
mtryTrials <- train(x_train_class, y_train_class, method = "rf", tuneLength = 6,
               trControl = cvControl, trace = FALSE)
mtryTrials
Random Forest 

537 samples
  7 predictor
  7 classes: 'A', 'B', 'C', 'D', 'E', 'F', 'G' 

No pre-processing
Resampling: Cross-Validated (10 fold) 
Summary of sample sizes: 483, 483, 485, 482, 483, 483, ... 
Resampling results across tuning parameters:

  mtry  Accuracy   Kappa    
  2     0.6554400  0.5830507
  3     0.6442214  0.5705064
  4     0.6423683  0.5690544
  5     0.6368451  0.5619084
  6     0.6443250  0.5713786
  7     0.6387008  0.5650792

Accuracy was used to select the optimal model using the largest value.
The final value used for the model was mtry = 2.
In [47]:
rf.c.opti <- randomForest::randomForest(x = x_train_class, y = y_train_class,
                         xtest = x_test_class, ytest = y_test_class,
                         mtry=mtryTrials$bestTune$mtry,importance = TRUE)

Résultats du modèle optimisé :

In [48]:
rf.c.opti.pred <- rf.c.opti$test$predicted
Tab_score_class <-Compute_Error(rf.c.opti.pred,quali=TRUE,name_model = "Random Forest Opti");Tab_score_class
Tab_F1_class <- Compute_F1(rf.c.opti.pred,quali=TRUE,name_model = "RF_Opti");
Table de contingence de Random Forest Opti : 
        observations
pred.reg  A  B  C  D  E  F  G
       A 55 21  3  0  0  0  0
       B  6 12  7  0  0  0  0
       C  0  4 10  0  0  0  0
       D  0  2  2 14  2  0  0
       E  0  0  0  6 15  7  0
       F  0  0  0  0 11 14  4
       G  0  0  0  0  0 12 24
ModeleAccuracyMAEMacroF1
Naive Logit 0.5194805 0.5714286 0.4594970
Penalised Logit 0.5021645 0.6363636 0.4465665
Tree opti 0.5887446 0.4545455 0.5638029
Random Forest 0.6233766 0.3982684 0.5929760
Random Forest Opti0.6233766 0.3982684 0.5933710

Commentaire :

Le modèle Random Forest retourne jusqu'ici les meilleurs résultats. La MAE est faible et les métriques F1 semblent être élevées.

On voit donc que l'apport du modèle optimisé pour le nombre de variables retenues aléatoirement à chaque noeud (mtry) n'est pas forcément significatif et ajoute un cout numérique supplémentaire.


e. Boosting

Mise en place du modèle :

In [49]:
defaultW <- getOption("warn")
options(warn = -1)

gbm.c <- gbm(Energy.efficiency ~ ., data=train_class,
        distribution = "multinomial",cv.folds = 10,shrinkage = .01,
       n.minobsinnode = 10, n.trees = 200, verbose=FALSE)

options(warn = defaultW)

Résultats :

In [50]:
p = predict.gbm(object = gbm.c,
                   newdata = test_class,
                   n.trees = 200,
                   type = "response")

gbm.c.pred = colnames(p)[apply(p, 1, which.max)]

Tab_score_class <- Compute_Error(gbm.c.pred,quali = TRUE, name_model = "GBM"); Tab_score_class
Tab_F1_class <- Compute_F1(gbm.c.pred,quali=TRUE,name_model = "GBM");
Table de contingence de GBM : 
        observations
pred.reg  A  B  C  D  E  F  G
       A 58 30  8  0  0  0  0
       B  1  2  2  0  0  0  0
       C  2  4  9  0  0  0  0
       D  0  3  2  8  2  1  0
       E  0  0  1 11 17 12  4
       F  0  0  0  1  6  6  0
       G  0  0  0  0  3 14 24
ModeleAccuracyMAEMacroF1
Naive Logit 0.5194805 0.5714286 0.4594970
Penalised Logit 0.5021645 0.6363636 0.4465665
Tree opti 0.5887446 0.4545455 0.5638029
Random Forest 0.6233766 0.3982684 0.5929760
Random Forest Opti0.6233766 0.3982684 0.5933710
GBM 0.5367965 0.5627706 0.4547098

Mise en place du modèle optimisé :

In [51]:
set.seed(2)
defaultW <- getOption("warn")
options(warn = -1)

cvControlRandom = trainControl(method = "cv", number = 10, search='random')
gbm.opti.c = train(x_train_class, y_train_class,method = "gbm", tuneLength = 8,
               trControl = cvControl, verbose=FALSE)

options(warn = defaultW)

Résultats :

In [52]:
gbm.opti.c.pred = predict(gbm.opti.c, newdata = x_test_class)
In [53]:
Tab_score_class <- Compute_Error(gbm.opti.c.pred,quali = TRUE, name_model = "GBM opti"); Tab_score_class
Tab_F1_class <- Compute_F1(gbm.opti.c.pred,quali=TRUE,name_model = "GBM opti")
Table de contingence de GBM opti : 
        observations
pred.reg  A  B  C  D  E  F  G
       A 53 19  3  0  0  0  0
       B  8 12  9  0  0  0  0
       C  0  5  7  0  1  0  0
       D  0  3  2 11  6  0  0
       E  0  0  0  8 10  6  0
       F  0  0  1  1  9 15  7
       G  0  0  0  0  2 12 21
ModeleAccuracyMAEMacroF1
Naive Logit 0.5194805 0.5714286 0.4594970
Penalised Logit 0.5021645 0.6363636 0.4465665
Tree opti 0.5887446 0.4545455 0.5638029
Random Forest 0.6233766 0.3982684 0.5929760
Random Forest Opti0.6233766 0.3982684 0.5933710
GBM 0.5367965 0.5627706 0.4547098
GBM opti 0.5584416 0.4935065 0.5088557

Commentaire :

Que ce soit le modèle optimisé ou non optimisé, utiliser la méthode de Gradient Boosting ne semble pas apporter de réelles améliorations par rapport à Random Forest, même avec des paramètres optimisés.


f. Support Vector Machine (SVM)

Mise en place du modèle :

In [54]:
svm.c=svm(Energy.efficiency~., data=train_class, kernel="radial")

Résultats :

In [55]:
svm.c.pred=predict(svm.c, newdata=test_class, type = "class")
Tab_score_class <- Compute_Error(svm.c.pred,quali = TRUE, name_model = "SVM"); Tab_score_class
Tab_F1_class <- Compute_F1(svm.c.pred,quali=TRUE,name_model = "SVM")
Table de contingence de SVM : 
        observations
pred.reg  A  B  C  D  E  F  G
       A 53 25 10  0  0  0  0
       B  6  6  4  0  0  0  0
       C  2  5  5  0  0  0  0
       D  0  3  3 12  7  5  3
       E  0  0  0  4  9  8  3
       F  0  0  0  2 10  9  3
       G  0  0  0  2  2 11 19
ModeleAccuracyMAEMacroF1
Naive Logit 0.5194805 0.5714286 0.4594970
Penalised Logit 0.5021645 0.6363636 0.4465665
Tree opti 0.5887446 0.4545455 0.5638029
Random Forest 0.6233766 0.3982684 0.5929760
Random Forest Opti0.6233766 0.3982684 0.5933710
GBM 0.5367965 0.5627706 0.4547098
GBM opti 0.5584416 0.4935065 0.5088557
SVM 0.4891775 0.6709957 0.4216265

Mise en place du modèle optimisé :

In [56]:
defaultW <- getOption("warn")
options(warn = -1)

svm.c.tune <- tune.svm(Energy.efficiency~., data=train_class, type = "C",
                       kernel=c('polynomial','radial','linear'), degree=c(2,3,4), 
                       cost=c(1,1.25,1.5,1.75,2),gamma = seq(0.1, 2, by = 0.2))

options(warn = defaultW)
In [57]:
(svm.c.tune$best.parameters)
degreegammacost
372 0.5 1.25

Résultats du modèle optimisé :

In [58]:
svm.c.opti=predict(svm.c.tune$best.model, newdata=test_class, type = "class")
In [59]:
Tab_score_class <- Compute_Error(svm.c.opti, quali = TRUE, name_model = "SVM Opti"); 
Tab_score_class
Tab_F1_class <- Compute_F1(svm.c.opti,quali=TRUE,name_model = "SVM Opti")
Table de contingence de SVM Opti : 
        observations
pred.reg  A  B  C  D  E  F  G
       A 50 23  7  0  0  0  0
       B  8  8  9  0  0  0  0
       C  3  5  4  1  0  0  0
       D  0  3  2 10  6  1  0
       E  0  0  0  5  9  6  3
       F  0  0  0  2 11 12  5
       G  0  0  0  2  2 14 20
ModeleAccuracyMAEMacroF1
Naive Logit 0.5194805 0.5714286 0.4594970
Penalised Logit 0.5021645 0.6363636 0.4465665
Tree opti 0.5887446 0.4545455 0.5638029
Random Forest 0.6233766 0.3982684 0.5929760
Random Forest Opti0.6233766 0.3982684 0.5933710
GBM 0.5367965 0.5627706 0.4547098
GBM opti 0.5584416 0.4935065 0.5088557
SVM 0.4891775 0.6709957 0.4216265
SVM Opti 0.4891775 0.6190476 0.4291337

Commentaire :

Tout comme GBM, la SVM n'apporte pas de meilleurs résultats de prédiction, la coût numérique de ces méthodes et de leur optimisation est d'ailleurs très élevé.


Bilan - Classification

In [60]:
cat("Bilan des scores de classification  : \n " )
Tab_score_class
Bilan des scores de classification  : 
 
ModeleAccuracyMAEMacroF1
Naive Logit 0.5194805 0.5714286 0.4594970
Penalised Logit 0.5021645 0.6363636 0.4465665
Tree opti 0.5887446 0.4545455 0.5638029
Random Forest 0.6233766 0.3982684 0.5929760
Random Forest Opti0.6233766 0.3982684 0.5933710
GBM 0.5367965 0.5627706 0.4547098
GBM opti 0.5584416 0.4935065 0.5088557
SVM 0.4891775 0.6709957 0.4216265
SVM Opti 0.4891775 0.6190476 0.4291337

Pour l'interprétation, se référer à la partie Conclusion.


2. Régression

Dans cette partie, nous allons traîter le problème comme un problème de regression. En effet, la variable Energy est une variable quantitative continue et c'est à partir de cette variable que les classes sont formées :

  • Les bâtiments ayant une consommation inférieure à 30 correspondent à une classe A
  • Les bâtiments ayant une consommation comprise entre 30 et 35 correspondent à une classe B
  • Les bâtiments ayant une consommation comprise entre 35 et 45 correspondent à une classe C
  • Les bâtiments ayant une consommation comprise entre 45 et 55 correspondent à une classe D
  • Les bâtiments ayant une consommation comprise entre 55 et 65 correspondent à une classe E
  • Les bâtiments ayant une consommation comprise entre 65 et 75 correspondent à une classe F
  • Les bâtiments ayant une consommation supérieure à 75 correspondent à une classe F

Ainsi nous allons prédire la variable continue Energy, puis attribuer une classe à la prédiction obtenue (selon la règle ci-dessus).

In [61]:
#New train and test sets : we remove Energy.efficiency
train_reg <- dplyr::select(train, -Energy.efficiency)
test_reg <- dplyr::select(test, -Energy.efficiency)
In [62]:
x_train_reg<-dplyr::select(train_reg, -Energy)
y_train_reg<-train_reg$Energy

x_test_reg<-dplyr::select(test_reg, -Energy)
y_test_reg<-test_reg$Energy
In [63]:
# permet de trouver la classe energetique à partir de la valeur numérique prédite

classify <- function(quantitative.predict){
    predict.class = (quantitative.predict)
    predict.class[quantitative.predict <=30] = 'A'
    predict.class[quantitative.predict >30  & quantitative.predict <= 35] = 'B'
    predict.class[quantitative.predict >35  & quantitative.predict <= 45] = 'C'
    predict.class[quantitative.predict >45  & quantitative.predict <= 55] = 'D'
    predict.class[quantitative.predict >55  & quantitative.predict <= 65] = 'E'
    predict.class[quantitative.predict >65  & quantitative.predict <= 75] = 'F'
    predict.class[quantitative.predict > 75] = 'G'   
    return(predict.class)
}

a. Regression linéaire

Mise en place du modèle :

In [64]:
glm.r <- glm(Energy ~., data=train_reg)

options(repr.plot.width = 8, repr.plot.height = 4) 
par(mfrow=c(1,2))
plot(glm.r, which=1:2)

Commentaires : Le graphe précédent nous permet de visualiser si les hypothèses de linéarité et d'homoscédasticité du modèle.

Les résidus sont centrés mais il semble y avoir une dépendance entre la valeur prédite et l'erreur (quand valeur prédite élevée : erreur forte); graphe en forme de banane.

In [65]:
summary(glm.r)
Call:
glm(formula = Energy ~ ., data = train_reg)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-19.8860   -3.9450   -0.3262    4.1105   25.5253  

Coefficients:
                             Estimate Std. Error t value Pr(>|t|)    
(Intercept)                 106.65152   35.73630   2.984  0.00297 ** 
Relative.compactness        -67.73059   20.83371  -3.251  0.00122 ** 
Surface.area                 -0.07005    0.03486  -2.009  0.04501 *  
Wall.area                     0.08454    0.01642   5.148 3.74e-07 ***
Overall.height.L             23.95255    1.98735  12.052  < 2e-16 ***
orientationNorth              0.30842    0.81037   0.381  0.70367    
orientationSouth              0.15928    0.82172   0.194  0.84638    
orientationWest               1.07865    0.84018   1.284  0.19977    
Glazing.area                 33.73610    2.44711  13.786  < 2e-16 ***
Glazing.area.distr55% North   4.30440    1.56811   2.745  0.00626 ** 
Glazing.area.distr55% East    4.94811    1.57624   3.139  0.00179 ** 
Glazing.area.distr55% South   3.82289    1.56007   2.450  0.01459 *  
Glazing.area.distr55% West    4.48659    1.56385   2.869  0.00429 ** 
Glazing.area.distrUniform     3.08239    1.56788   1.966  0.04983 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 45.67094)

    Null deviance: 206282  on 536  degrees of freedom
Residual deviance:  23886  on 523  degrees of freedom
AIC: 3591.9

Number of Fisher Scoring iterations: 2

Résultats :

In [66]:
glm.r.pred <- predict(glm.r, newdata=test_reg)
In [67]:
Tab_score_reg <- Compute_Error(glm.r.pred,quali = FALSE, name_model = "Regression Linéaire Naive");Tab_score_reg
Tab_F1_reg <- Compute_F1(glm.r.pred,quali=FALSE,name_model = "Reg Naive")
Table de contingence de Regression Linéaire Naive : 
        observations
pred.reg  A  B  C  D  E  F  G
       A 43 12  3  0  0  0  0
       B 13 17  9  0  0  0  0
       C  5  7  7  0  0  0  0
       D  0  3  3  4  2  0  0
       E  0  0  0 12 17 17  3
       F  0  0  0  4  9 15 19
       G  0  0  0  0  0  1  6
ModeleAccuracyMAEMacroF1RMSE
Regression Linéaire Naive0.4718615 0.6060606 0.4156379 7.292765

Commentaire :

Une régression linéaire naive a été réalisée : tous les paramètres ont été ajoutés au modèle.

La RMSE reste élevé pour distinguer toutes les classes. En effet, l'étendu des différentes classes est souvent égale à 10 (même parfois à 5 pour la classe B). La RMSE reste donc trop élevé lors de la séparation entre les différentes classes.

Les résultats obtenues avec la regression ne semble pas apporter une amélioration au modèle.


b. Régression linéaire pénalisée

Mise en place du modèle :

In [68]:
x.mat.train <- model.matrix(Energy ~ . - 1, data = train_reg) #permet de gérer les variables catégorielles
lasso.r <- glmnet(y = y_train_reg, x = x.mat.train)

options(repr.plot.width = 10, repr.plot.height = 7)
plot(lasso.r, xvar = "lambda", label = TRUE)
grid()
legend("bottomright", 
       legend = paste(1:ncol(x.mat.train), " - ", colnames(x.mat.train)), ncol=2, cex=0.75)

# choix du paramètre de régularisation par validation croisée

lasso.r.cv <- cv.glmnet(y = y_train_reg, x = x.mat.train)
plot(lasso.r.cv)
In [69]:
#Coefficients du modèle retenus par la pénalisation LASSO
coef(lasso.r.cv, s = "lambda.1se")
15 x 1 sparse Matrix of class "dgCMatrix"
                                        1
(Intercept)                  2.992255e+01
Relative.compactness         .           
Surface.area                 .           
Wall.area                    7.957880e-02
Overall.height3.5           -3.091465e+01
Overall.height7              1.485976e-11
orientationNorth             .           
orientationSouth             .           
orientationWest              .           
Glazing.area                 3.042508e+01
Glazing.area.distr55% North  .           
Glazing.area.distr55% East   .           
Glazing.area.distr55% South  .           
Glazing.area.distr55% West   .           
Glazing.area.distrUniform    .           

Résultats :

In [70]:
x.mat.test <- model.matrix(Energy ~ . - 1, data = test_reg)
lasso.r.cv.pred <- predict(lasso.r.cv, s = "lambda.1se", newx = x.mat.test) #valeur predite sur le jeu de train
In [71]:
Tab_score_reg <- Compute_Error(lasso.r.cv.pred,quali = FALSE, name_model = "Lasso penalisation");Tab_score_reg
Tab_F1_reg <- Compute_F1(lasso.r.cv.pred,quali=FALSE,name_model = "Reg LASSO")
Table de contingence de Lasso penalisation : 
        observations
pred.reg  A  B  C  D  E  F  G
       A 40 11  2  0  0  0  0
       B 15 15 10  0  0  0  0
       C  6 10  7  0  0  0  0
       D  0  1  1  1  0  0  0
       E  0  2  2 17 18 16  5
       F  0  0  0  2 10 17 17
       G  0  0  0  0  0  0  6
ModeleAccuracyMAEMacroF1RMSE
Regression Linéaire Naive0.4718615 0.6060606 0.4156379 7.292765
Lasso penalisation 0.4502165 0.6450216 0.3817115 7.823848

Rappel : La pénalisation LASSO permet de diminuer le nombre de variable utilisé tout en dégradant très peu l'information. Dans ce modèle une contrainte est ajouté sur les paramètres. La pénalisation LASSO a donc une RMSE supérieur au modèle linéaire car le modèle linéaire minimise uniquement les moindres carrés. Nous avons donc un moins bon modèle mais moins de variables sont nécessaires pour le construire.

c. Régression linéaire avec sélection de variable par le critère Akaike Information Criterion (AIC) :

Mise en place du modèle :

In [72]:
#Test sur le modèle sans les intéractions
glm.r <- glm(Energy~ ., data=train_reg)
glm.r.step <- step(glm.r, direction = "backward", trace=0)
anova(glm.r.step, test='F')
DfDevianceResid. DfResid. DevFPr(>F)
NULLNA NA 536 206282.13 NA NA
Relative.compactness 1 79795.9899 535 126486.14 1750.868619 1.727430e-169
Surface.area 1 15724.8126 534 110761.32 345.030883 1.353660e-59
Wall.area 1 66075.8896 533 44685.43 1449.824755 2.789812e-153
Overall.height 1 7727.8702 532 36957.56 169.563475 8.509584e-34
Glazing.area 1 12440.5465 531 24517.02 272.968133 1.063854e-49
Glazing.area.distr 5 544.5219 526 23972.50 2.389559 3.693511e-02

Résultats :

In [73]:
glm.r.step.pred <- predict(glm.r.step, newdata=test_reg)
Tab_score_reg <- Compute_Error(glm.r.step.pred,quali = FALSE, name_model = "Regression Linéaire Selection Variable");Tab_score_reg
Tab_F1_reg <- Compute_F1(glm.r.step.pred,quali=FALSE,name_model = "Reg AIC")
Table de contingence de Regression Linéaire Selection Variable : 
        observations
pred.reg  A  B  C  D  E  F  G
       A 44 12  4  0  0  0  0
       B 12 16  8  0  0  0  0
       C  5  8  7  0  0  0  0
       D  0  3  3  5  2  0  0
       E  0  0  0 10 17 17  3
       F  0  0  0  5  9 14 19
       G  0  0  0  0  0  2  6
ModeleAccuracyMAEMacroF1RMSE
Regression Linéaire Naive 0.4718615 0.6060606 0.4156379 7.292765
Lasso penalisation 0.4502165 0.6450216 0.3817115 7.823848
Regression Linéaire Selection Variable0.4718615 0.6147186 0.4181385 7.278510

Mise en place du modèle quadratique :

In [74]:
glm.2.r <- glm(Energy~ .^2, data=train_reg)
# du critère d'Akaïke par méthode descendante
glm.2.r.step <- step(glm.2.r, direction = "backward", trace=0)
In [75]:
#On réalise un test de fischer sur le résultat pour confirmer que toutes les variables retenues sont bien significatives
anova(glm.2.r.step, test='F')
DfDevianceResid. DfResid. DevFPr(>F)
NULLNA NA 536 206282.13 NA NA
Relative.compactness 1 79795.98993 535 126486.14 2084.9638165 3.177342e-181
Surface.area 1 15724.81256 534 110761.32 410.8685816 2.974337e-67
Wall.area 1 66075.88960 533 44685.43 1726.4757175 6.640534e-165
Overall.height 1 7727.87016 532 36957.56 201.9190397 8.839812e-39
Glazing.area 1 12440.54646 531 24517.02 325.0550466 1.914205e-56
Glazing.area.distr 5 544.52186 526 23972.50 2.8455274 1.515766e-02
Relative.compactness:Surface.area 1 722.54386 525 23249.95 18.8791168 1.685023e-05
Relative.compactness:Wall.area 1 403.61526 524 22846.34 10.5459336 1.242009e-03
Relative.compactness:Overall.height 1 734.07126 523 22112.27 19.1803124 1.447307e-05
Relative.compactness:Glazing.area 1 1591.48278 522 20520.78 41.5833428 2.649308e-10
Surface.area:Overall.height 1 76.38985 521 20444.39 1.9959660 1.583357e-01
Surface.area:Glazing.area.distr 5 108.52184 516 20335.87 0.5671065 7.252539e-01
Wall.area:Overall.height 1 79.47461 515 20256.40 2.0765666 1.501973e-01
Wall.area:Glazing.area 1 243.72448 514 20012.67 6.3681986 1.192468e-02
Wall.area:Glazing.area.distr 5 111.34047 509 19901.33 0.5818359 7.139460e-01
Overall.height:Glazing.area.distr 5 612.18187 504 19289.15 3.1991007 7.458697e-03

Résultats :

In [76]:
glm.2.r.step.pred <- predict(glm.2.r.step, newdata=test_reg)
In [77]:
Tab_score_reg <- Compute_Error(glm.2.r.step.pred,quali = FALSE, name_model = "Regression Linéaire Quadratique Selection Variable");Tab_score_reg
Tab_F1_reg <- Compute_F1(glm.2.r.step.pred,quali=FALSE,name_model = "Reg AIC quadra")
Table de contingence de Regression Linéaire Quadratique Selection Variable : 
        observations
pred.reg  A  B  C  D  E  F  G
       A 44 14  3  0  0  0  0
       B 15 19 14  0  0  0  0
       C  2  4  3  0  0  0  0
       D  0  0  0  9  5  0  0
       E  0  2  2  7 13 10  3
       F  0  0  0  4 10 19 13
       G  0  0  0  0  0  4 12
ModeleAccuracyMAEMacroF1RMSE
Regression Linéaire Naive 0.4718615 0.6060606 0.4156379 7.292765
Lasso penalisation 0.4502165 0.6450216 0.3817115 7.823848
Regression Linéaire Selection Variable 0.4718615 0.6147186 0.4181385 7.278510
Regression Linéaire Quadratique Selection Variable0.5151515 0.5627706 0.4725029 7.148471

Commentaires : Deux regressions linéaires avec selection de variables ont été réalisées : le modèle initial et un modèle quadratique. La selection de variables (ici méthode descendante) permet de garder les variables les plus intéressantes au sens du critère AIC.

En comparant les différents scores, la regression quadratique avec la selection de variables est très légèrement meilleure.


d. Arbre de décision optimal

Mise en place du modèle :

In [78]:
options(repr.plot.width = 15, repr.plot.height = 10)
tree.r=rpart(Energy~.,data=train_reg,control=rpart.control(cp=0.001))
In [79]:
xmat=xpred.rpart(tree.r, xval = 10)
#levels(datappr[,"Survived"]) <- c(1,2)
xerr <- (xmat-as.numeric(y_train_reg))^2
CVerr=apply(xerr,2,sum)
#CVerr  #    CP           erreur
cp_opti = as.numeric(attributes(which.min(CVerr))$names); 

cat('cp opti :', cp_opti)

tree.opti.r=rpart(Energy~.,data=train_reg,control=rpart.control(cp=cp_opti))
options(repr.plot.width = 15, repr.plot.height = 10)
plot(tree.opti.r)
text(tree.opti.r)
cp opti : 0.001072396

Commentaire : On remarque que, comme pour le modèle de classification, la variable Overall.height a un grand impact sur le modèle.

Résultats :

In [80]:
tree.opti.r.pred <- predict(tree.opti.r, newdata = test_reg)
Tab_score_reg <- Compute_Error(tree.opti.r.pred,quali = FALSE, name_model = "Tree Opti");Tab_score_reg
Tab_F1_reg <- Compute_F1(tree.opti.r.pred,quali=FALSE,name_model = "TreeOpti Reg")
Table de contingence de Tree Opti : 
        observations
pred.reg  A  B  C  D  E  F  G
       A 46  8  2  0  0  0  0
       B 12 20  4  0  0  0  0
       C  3 11 16  1  0  0  0
       D  0  0  0 14  4  0  0
       E  0  0  0  5 11  5  0
       F  0  0  0  0 13 19  7
       G  0  0  0  0  0  9 21
ModeleAccuracyMAEMacroF1RMSE
Regression Linéaire Naive 0.4718615 0.6060606 0.4156379 7.292765
Lasso penalisation 0.4502165 0.6450216 0.3817115 7.823848
Regression Linéaire Selection Variable 0.4718615 0.6147186 0.4181385 7.278510
Regression Linéaire Quadratique Selection Variable0.5151515 0.5627706 0.4725029 7.148471
Tree Opti 0.6363636 0.3852814 0.6230242 5.458206

Commentaire :

Tout comme précedemment en régression, nous avons une amélioration significative des résultats avec un arbre de décision (méthode non linéaire).


e. Random Forest

Mise en place du modèle :

In [81]:
set.seed(69)
rf.r <-randomForest::randomForest(x = x_train_reg, y = y_train_reg,
                         xtest = x_test_reg, ytest = y_test_reg,
                         importance = TRUE)

Résultats :

In [82]:
rf.r.pred <- rf.r$test$predicted
Tab_score_reg <- Compute_Error(rf.r.pred,quali = FALSE, name_model = "Random Forest");Tab_score_reg
Tab_F1_reg <- Compute_F1(rf.r.pred,quali=FALSE,name_model = "Reg RF")
Table de contingence de Random Forest : 
        observations
pred.reg  A  B  C  D  E  F  G
       A 47 11  1  0  0  0  0
       B 13 24  7  0  0  0  0
       C  1  2 11  0  0  0  0
       D  0  2  3 11  1  0  0
       E  0  0  0  9 18  8  1
       F  0  0  0  0  9 18  9
       G  0  0  0  0  0  7 18
ModeleAccuracyMAEMacroF1RMSE
Regression Linéaire Naive 0.4718615 0.6060606 0.4156379 7.292765
Lasso penalisation 0.4502165 0.6450216 0.3817115 7.823848
Regression Linéaire Selection Variable 0.4718615 0.6147186 0.4181385 7.278510
Regression Linéaire Quadratique Selection Variable0.5151515 0.5627706 0.4725029 7.148471
Tree Opti 0.6363636 0.3852814 0.6230242 5.458206
Random Forest 0.6363636 0.3852814 0.6186910 5.155054

Commentaire :

Les résultats obtenus avec Random Forest sont similaires à l'arbre de décision optimisé. Cependant, la RMSE est meilleur pour Random Forest.

Tout comme précemment, nous allons tenter d'améliorer le modèle Random Forest en optimisant le paramètre mtry correspondant au nombre de variables testées à chaque division.

Mise en place du modèle optimisé :

In [83]:
set.seed(69)
cvControl <- trainControl(method = "cv", number = 10)
mtryTrials <- train(train_reg, y_train_reg, method = "rf", tuneLength = 7,
               trControl = cvControl, trace = FALSE)
In [84]:
set.seed(69)
defaultW <- getOption("warn")
options(warn = -1)
rf.opti.r<-randomForest::randomForest(x = x_train_reg, y = y_train_reg,
                         xtest = x_test_reg, ytest = y_test_reg,
                         mtry=mtryTrials$bestTune$mtry,importance = TRUE)

options(warn = defaultW)

Résultats :

In [85]:
rf.opti.r.pred <- rf.opti.r$test$predicted
Tab_score_reg <- Compute_Error(rf.opti.r.pred,quali = FALSE, name_model = "Random Forest Opti");Tab_score_reg
Tab_F1_reg <- Compute_F1(rf.opti.r.pred,quali=FALSE,name_model = "Reg RF Opti")
Table de contingence de Random Forest Opti : 
        observations
pred.reg  A  B  C  D  E  F  G
       A 46  8  1  0  0  0  0
       B 14 26  8  0  0  0  0
       C  1  3 11  0  0  0  0
       D  0  2  2 13  2  0  0
       E  0  0  0  7 15  7  0
       F  0  0  0  0 11 18 11
       G  0  0  0  0  0  8 17
ModeleAccuracyMAEMacroF1RMSE
Regression Linéaire Naive 0.4718615 0.6060606 0.4156379 7.292765
Lasso penalisation 0.4502165 0.6450216 0.3817115 7.823848
Regression Linéaire Selection Variable 0.4718615 0.6147186 0.4181385 7.278510
Regression Linéaire Quadratique Selection Variable0.5151515 0.5627706 0.4725029 7.148471
Tree Opti 0.6363636 0.3852814 0.6230242 5.458206
Random Forest 0.6363636 0.3852814 0.6186910 5.155054
Random Forest Opti 0.6320346 0.3852814 0.6161488 5.189220

Commentaire :

Le modèle Random Forest retourne jusqu'ici les meilleurs résultats en terme de RMSE.

Cependant, par rapport aux autres métriques, les modèles Random Forest et d'arbre de décision sont équivalents.


f. Boosting

Mise en place du modèle :

In [86]:
gbm.r=gbm(Energy~., data=train_reg,distribution="gaussian",n.trees=2000, cv.folds=10,
        n.minobsinnode = 5,shrinkage=0.05,verbose=FALSE)
plot(gbm.r$cv.error)
In [87]:
# nombre optimal d'itérations par valiation croisée
best.iter=gbm.perf(gbm.r,method="cv")
In [88]:
gbm.r.pred=predict(gbm.r,newdata=test_reg,n.trees=best.iter)

Tab_score_reg <- Compute_Error(gbm.r.pred,quali = FALSE, name_model = "GBM Naive");Tab_score_reg
Tab_F1_reg <- Compute_F1(gbm.r.pred,quali=FALSE,name_model = "Reg GBM")
Table de contingence de GBM Naive : 
        observations
pred.reg  A  B  C  D  E  F  G
       A 45  9  2  0  0  0  0
       B 12 18  7  0  0  0  0
       C  4 10 12  0  0  0  0
       D  0  2  1  8  4  0  0
       E  0  0  0 12 17  9  1
       F  0  0  0  0  7 19  9
       G  0  0  0  0  0  5 18
ModeleAccuracyMAEMacroF1RMSE
Regression Linéaire Naive 0.4718615 0.6060606 0.4156379 7.292765
Lasso penalisation 0.4502165 0.6450216 0.3817115 7.823848
Regression Linéaire Selection Variable 0.4718615 0.6147186 0.4181385 7.278510
Regression Linéaire Quadratique Selection Variable0.5151515 0.5627706 0.4725029 7.148471
Tree Opti 0.6363636 0.3852814 0.6230242 5.458206
Random Forest 0.6363636 0.3852814 0.6186910 5.155054
Random Forest Opti 0.6320346 0.3852814 0.6161488 5.189220
GBM Naive 0.5930736 0.4458874 0.5674609 5.529532

Mise en place du modèle optimisé:

In [89]:
cvControlRandom = trainControl(method = "cv", number = 10, search='random')
gbm.tune.r = train(x_train_reg, y_train_reg, method='gbm', tuneLength=4, trControl=cvControlRandom, verbose=FALSE)

Résultats :

In [90]:
gbm.tune.r.pred = predict(gbm.tune.r, newdata = x_test_reg)
In [91]:
Tab_score_reg <- Compute_Error(gbm.tune.r.pred, quali = FALSE, 
                               name_model = "GBM Opti");Tab_score_reg
Tab_F1_reg <- Compute_F1(gbm.tune.r.pred,quali=FALSE,name_model = "Reg GBM Opti")
Table de contingence de GBM Opti : 
        observations
pred.reg  A  B  C  D  E  F  G
       A 47 17  3  0  0  0  0
       B 10 14  5  0  0  0  0
       C  3  8 12  1  0  0  0
       D  1  0  2  9  3  0  0
       E  0  0  0 10 20  6  2
       F  0  0  0  0  5 19  5
       G  0  0  0  0  0  8 21
ModeleAccuracyMAEMacroF1RMSE
Regression Linéaire Naive 0.4718615 0.6060606 0.4156379 7.292765
Lasso penalisation 0.4502165 0.6450216 0.3817115 7.823848
Regression Linéaire Selection Variable 0.4718615 0.6147186 0.4181385 7.278510
Regression Linéaire Quadratique Selection Variable0.5151515 0.5627706 0.4725029 7.148471
Tree Opti 0.6363636 0.3852814 0.6230242 5.458206
Random Forest 0.6363636 0.3852814 0.6186910 5.155054
Random Forest Opti 0.6320346 0.3852814 0.6161488 5.189220
GBM Naive 0.5930736 0.4458874 0.5674609 5.529532
GBM Opti 0.6147186 0.4285714 0.5911386 5.859907

Commentaire :

Tout comme avec la classification, la méthode de Gradient Boosting ne semble pas apporter de réelles améliorations par rapport à Random Forest, même avec des paramètres optimisés.


g. SVM

Mise en place du modèle optimisé:

In [92]:
svm.tune.r = tune.svm(Energy~.,data=train_reg,
                       cost=c(1,1.25,1.5,1.75,2),gamma = seq(0.1, 2, by = 0.2))

Commentaires :

Etant limité par les performances de nos machines, nous ne pouvons pas tester plusieurs noyaux pour optimiser le modèle SVM en régression.

Résultats :

In [93]:
svm.opti.r.pred=predict(svm.tune.r$best.model,newdata=test_reg)
Tab_score_reg <- Compute_Error(svm.opti.r.pred,quali = FALSE, name_model = "SVM Opti");Tab_score_reg
Tab_F1_reg <- Compute_F1(svm.opti.r.pred,quali=FALSE,name_model = "Reg SVM Opti")
Table de contingence de SVM Opti : 
        observations
pred.reg  A  B  C  D  E  F  G
       A 39 12  1  0  0  0  0
       B 18 16 10  0  0  0  0
       C  4  9  8  1  0  0  0
       D  0  1  3  8  4  0  0
       E  0  1  0  8 12  6  0
       F  0  0  0  3 12 18 14
       G  0  0  0  0  0  9 14
ModeleAccuracyMAEMacroF1RMSE
Regression Linéaire Naive 0.4718615 0.6060606 0.4156379 7.292765
Lasso penalisation 0.4502165 0.6450216 0.3817115 7.823848
Regression Linéaire Selection Variable 0.4718615 0.6147186 0.4181385 7.278510
Regression Linéaire Quadratique Selection Variable0.5151515 0.5627706 0.4725029 7.148471
Tree Opti 0.6363636 0.3852814 0.6230242 5.458206
Random Forest 0.6363636 0.3852814 0.6186910 5.155054
Random Forest Opti 0.6320346 0.3852814 0.6161488 5.189220
GBM Naive 0.5930736 0.4458874 0.5674609 5.529532
GBM Opti 0.6147186 0.4285714 0.5911386 5.859907
SVM Opti 0.4978355 0.5497835 0.4741817 6.745372

Commentaire :

Nous n'obtenons pas de meilleur résultat avec SVM, de plus cette méthode est très couteuse.


Bilan - Régression

In [94]:
cat("Bilan des scores en passant par la Regression  : \n " )
Tab_score_reg
Bilan des scores en passant par la Regression  : 
 
ModeleAccuracyMAEMacroF1RMSE
Regression Linéaire Naive 0.4718615 0.6060606 0.4156379 7.292765
Lasso penalisation 0.4502165 0.6450216 0.3817115 7.823848
Regression Linéaire Selection Variable 0.4718615 0.6147186 0.4181385 7.278510
Regression Linéaire Quadratique Selection Variable0.5151515 0.5627706 0.4725029 7.148471
Tree Opti 0.6363636 0.3852814 0.6230242 5.458206
Random Forest 0.6363636 0.3852814 0.6186910 5.155054
Random Forest Opti 0.6320346 0.3852814 0.6161488 5.189220
GBM Naive 0.5930736 0.4458874 0.5674609 5.529532
GBM Opti 0.6147186 0.4285714 0.5911386 5.859907
SVM Opti 0.4978355 0.5497835 0.4741817 6.745372

Pour l'interprétation, se référer à la partie Conclusion.


D. Conclusion

In [95]:
cat("Bilan des scores - Classification : \n " )
Tab_score_class
Bilan des scores - Classification : 
 
ModeleAccuracyMAEMacroF1
Naive Logit 0.5194805 0.5714286 0.4594970
Penalised Logit 0.5021645 0.6363636 0.4465665
Tree opti 0.5887446 0.4545455 0.5638029
Random Forest 0.6233766 0.3982684 0.5929760
Random Forest Opti0.6233766 0.3982684 0.5933710
GBM 0.5367965 0.5627706 0.4547098
GBM opti 0.5584416 0.4935065 0.5088557
SVM 0.4891775 0.6709957 0.4216265
SVM Opti 0.4891775 0.6190476 0.4291337
In [96]:
cat("Bilan des scores - Regression  : \n " )
Tab_score_reg
Bilan des scores - Regression  : 
 
ModeleAccuracyMAEMacroF1RMSE
Regression Linéaire Naive 0.4718615 0.6060606 0.4156379 7.292765
Lasso penalisation 0.4502165 0.6450216 0.3817115 7.823848
Regression Linéaire Selection Variable 0.4718615 0.6147186 0.4181385 7.278510
Regression Linéaire Quadratique Selection Variable0.5151515 0.5627706 0.4725029 7.148471
Tree Opti 0.6363636 0.3852814 0.6230242 5.458206
Random Forest 0.6363636 0.3852814 0.6186910 5.155054
Random Forest Opti 0.6320346 0.3852814 0.6161488 5.189220
GBM Naive 0.5930736 0.4458874 0.5674609 5.529532
GBM Opti 0.6147186 0.4285714 0.5911386 5.859907
SVM Opti 0.4978355 0.5497835 0.4741817 6.745372

Pour rappel, nous pouvons considérer un modèle comme étant performant si il a une forte accuracy, un faible MAE et un MacroF1 élevé. Dans le cas d'une régression, on considère la RMSE comme étant un bon indicateur.

En classification, le meilleur modèle obtenu est celui de Random Forest : il est gagnant sur toutes les métriques.

En régression, nous observons une équivalence des modèles Random Forest et Arbres de décision : Random Forest a une RMSE plus faible mais une Macro F1 moins bonne que celle de l'arbre de décision.

A première vue l'accuracy ne semble pas élevé, mais nous pouvons nuancer cette constatation par le fait que la classification se fait selon 7 classes. Aléatoirement, nous avons 1 chance sur 7 de choisir la bonne classe soit une "accuracy" aléatoire de 0.15. Ici notre modèle est donc nettement supérieure à cette probabilité.

De plus, les erreurs de classification sont souvent réalisées dans les classes voisines (cf MAE).

En prenant en compte tous ces éléments nos modèles semblent relativement pertinents.


Analyse des scores F1

Les graphiques suivants exposent les scores F1 de chaque classe énergétique pour chaque modèle implémenté.

In [97]:
l = nrow(Tab_F1_class)
c = ncol(Tab_F1_class)
gras = 3
plot(Tab_F1_class$Modele,rep(-1,l),type='l', ylim=c(0,1),
     xlab="Modele", ylab= "F1 Score", main = "Score F1 pour les différentes classes énergétiques pour les modèles de classification")
grid()
points(Tab_F1_class$Modele,Tab_F1_class$A,  col=1, lwd=gras)
points(Tab_F1_class$Modele,Tab_F1_class$B, col=2, lwd = gras)
points(Tab_F1_class$Modele,Tab_F1_class$C, col=3, lwd = gras)
points(Tab_F1_class$Modele,Tab_F1_class$D, col=4, lwd = gras)
points(Tab_F1_class$Modele,Tab_F1_class$E, col=5, lwd = gras)
points(Tab_F1_class$Modele,Tab_F1_class$F, col=6, lwd = gras)
points(Tab_F1_class$Modele,Tab_F1_class$G, col=7, lwd = gras)
legend(x='bottomright', legend=c("A","B","C","D","E","F","G"), pch='o', col=1:(c-1))
In [98]:
l = nrow(Tab_F1_reg)
c = ncol(Tab_F1_reg)
gras = 3
plot(Tab_F1_reg$Modele,rep(-1,l),type='l', ylim=c(0,1),
     xlab="Modele", ylab= "F1 Score", main = "Score F1 pour les différentes classes énergétiques pour les modèles de régression")
grid()
points(Tab_F1_reg$Modele,Tab_F1_reg$A,  col=1, lwd=gras)
points(Tab_F1_reg$Modele,Tab_F1_reg$B, col=2, lwd = gras)
points(Tab_F1_reg$Modele,Tab_F1_reg$C, col=3, lwd = gras)
points(Tab_F1_reg$Modele,Tab_F1_reg$D, col=4, lwd = gras)
points(Tab_F1_reg$Modele,Tab_F1_reg$E, col=5, lwd = gras)
points(Tab_F1_reg$Modele,Tab_F1_reg$F, col=6, lwd = gras)
points(Tab_F1_reg$Modele,Tab_F1_reg$G, col=7, lwd = gras)
legend(x='bottomright', legend=c("A","B","C","D","E","F","G"), pch='o', col=1:(c-1))

En classification, ce sont les classes A et G qui sont les mieux prédites pour tous les modèles. Ce résultat est cohérent car ce sont des classes extrêmes avec de larges intervalles (A $\in [0,30[$ et G>75). En régression, cette constatation est également vraie hormis pour les modèles linéaires.

En classification, la classe B a du mal à être prédite. Ce résultat est aussi chohérent car l'intervalle de la classe B est inférieur aux autres (B $\in [30,35[$). Le risque de se tromper est donc augmenté. En régression, on ne distingue pas de classe majoritairement discriminée.

En mettant en parallèle les deux graphiques : il est interessant de constater pour Random Forest que les valeurs de macroF1 sont plus élevées en regression que en classification. Par ailleurs le delta en regression est moins important que celui en classification. On peut donc considérer le modèle de regression comme meilleur.

Une piste d'ouverture de notre étude consisterait à implémenter un algorithme permettant d'abord de prédire le type de consommation (Faible ou Forte) correspondant aux classes (A,B,C ou D,E,F,G) et ensuite d'appliquer un algo de classification sur ces données.


Lien entre les modèles et l'analyse statistique initiale :

Les résultats obtenues avec les différents modèles sont cohérents avec notre analyse exploratoire des données.

  • Répartition des classes

On voit tout d'abord qus les classes B,C,D sont très peu représentées dans le jeu de données initial et que la classe A domine. Cela a un impact direct sur la classification, les individus auront tendance à être prédits dans la classe A, même si ils sont dans la classe B ou la classe C. Nos modèles apprennent sur moins de données provenant de ses classes, tout cela combiné à leur situation de classes intermédiaires (donc plus difficiles à prédire car non extrêmes) entraine une réelle difficulté à vraiment segmenter les classes B et C, comme le montrent nos graphiques représentant les métriques F1 selon les classes.

Une solution pourrait être de choisir des jeux d'apprentissages avec plus d'individus de ces classes mals représentées.

  • Variables retenues dans les modèles

Ensuite on voit que dans nos modèles pénalisés les variables retenues sont, dans le cas de LASSO avec lambda=1SE, Overall.height, Glazing Area et Wall Area. Cela s'explique car les variables Wall.Area et Glazing.Area sont les variables les moins corrélées avec les autres dans notre matrice des corrélations. On ne peut donc pas s'en passer car leur information n'est liée à aucune autre variable et ne peut donc se retrouver ailleurs.

On observe aussi que la première variable servant à départager le premier noeud des arbres est overall.height. Dès les premiers graphiques montrant la hauteur selon la classe énergétique, on observait déjà que la variable Overall.height jouait un rôle important. En effet, les classes A,B et C sont toutes de hauteur 3,5m contre 7m pour les classes D,E,F et G. La deuxième variable intervenant dans les arbres pour la séparation est Glazing Area, variable portant à elle seule l'axe 2 de l'ACP (et expliquant 25% de la variance). Cette variable a donc évidemment un rôle important dans les modèles.

Enfin, le clustering ne mettait pas en évidence 7 clusters, ce qui laissait présager que la classification ne serait pas évidente. On remarque par ailleurs que la méthode du coude nous suggère 2 clusters (correspondant aux faibles et fortes consommations) et c'est réellement sur ce point que nous sommes bons : prédire une forte ou faible consommation (cf partie Bonus).

Enfin, les méthodes d'apprentissage non linéaires sont plus efficaces ici, même en regression.

Bonus - classification hierarchique

Nous avons vu de part le clustering et l'ACP, deux classes se distinguaient bien (ABC vs DEFG) nous allons donc mettre en place une classification hiérarchique :

1. Mise en place d'un modèle de prédiction binaire : Bonne ou mauvaise consommation ?

In [99]:
#On fait une nouvelle data_frame relabellisée : 
x_bin = x_fin
In [100]:
Conso = ifelse(x_fin$Energy<=45,"Faible","Forte") 
x_bin = cbind(x_bin,Conso)
x_bin$Conso <- factor(x_bin$Conso)
In [101]:
head(x_bin)
Relative.compactnessSurface.areaWall.areaOverall.heightorientationGlazing.areaGlazing.area.distrEnergyEnergy.efficiencyConso
0.9829276 530.4900 306.4846 7 North 1.609490e-02No Glazing 34.26394 B Faible
0.9835473 519.8724 299.7763 7 East 0.000000e+00No Glazing 34.58975 B Faible
0.9794535 516.1912 303.3744 7 South 0.000000e+00No Glazing 38.77805 C Faible
0.9777325 518.9241 292.8122 7 West 9.554434e-06No Glazing 37.94781 C Faible
0.9030294 552.9689 316.2361 7 North 0.000000e+00No Glazing 47.67586 D Forte
0.8909102 558.6037 314.9162 7 East 0.000000e+00No Glazing 41.90847 C Faible
In [102]:
train_bin <- x_bin[train_ind, ]
test_bin <- x_bin[-train_ind, ]


train_bin <- dplyr::select(train_bin, Relative.compactness,Surface.area, Wall.area,Overall.height,orientation, Glazing.area,Glazing.area.distr,Conso)
test_bin <- dplyr::select(test_bin, Relative.compactness, Surface.area,Wall.area,Overall.height,orientation, Glazing.area,Glazing.area.distr,Conso)
In [103]:
prems <- glm(Conso ~., data=train_bin, family = binomial(link = "logit"))
In [104]:
res.fit <- predict(prems, newdata = test_bin, type="response")
In [105]:
t1=table(res.fit>0.5,test_bin$Conso);t1
score = sum(diag(t1)/sum(t1))
print(score)
       
        Faible Forte
  FALSE    116     1
  TRUE       6   108
[1] 0.969697

2. Mise en place des modèles de prédiction pour les consos fortes et faibles

Séparation des données

In [106]:
train_bin <- x_bin[train_ind, ]
train_bin <- dplyr::select(train_bin,-Energy.efficiency)

test_bin <- x_bin[-train_ind, ]
test_bin <- dplyr::select(test_bin,-Energy.efficiency)

#juste : on sépare les forts et les faibles reels pour entrainer nos modèles
train_bin_Forte <-train_bin[which(train_bin$Conso=='Forte'),]
train_bin_Faible <-train_bin[which(train_bin$Conso=='Faible'),]
#vraies valeurs correspondant à nos predictions de bonne ou mauvaise classe
test_bin_Forte <- test_bin[which(res.fit > 0.5),]
test_bin_Faible <- test_bin[which(res.fit < 0.5),]

#On enlève la variable Conso : 
train_bin_Forte <- dplyr::select(train_bin_Forte,-Conso)
train_bin_Faible <- dplyr::select(train_bin_Faible,-Conso)

test_bin_Forte <- dplyr::select(test_bin_Forte, -Conso)
test_bin_Faible <- dplyr::select(test_bin_Faible,-Conso)


x_train_bin_Forte <- dplyr::select(train_bin_Forte, -Energy)
y_train_bin_Forte <- train_bin_Forte$Energy

x_test_bin_Forte <- dplyr::select(test_bin_Forte, -Energy)
y_test_bin_Forte <- test_bin_Forte$Energy

x_train_bin_Faible <- dplyr::select(train_bin_Faible, -Energy)
y_train_bin_Faible <- train_bin_Faible$Energy

x_test_bin_Faible <- dplyr::select(test_bin_Faible, -Energy)
y_test_bin_Faible <- test_bin_Faible$Energy

Modele de prediction : fortes consommations

In [107]:
rf.h.forte <- randomForest(x = x_train_bin_Forte, y = y_train_bin_Forte ,
                      xtest = x_test_bin_Forte, ytest = y_test_bin_Forte,
                     ntree = 500,do.trace = 50, importance = TRUE)

rf.h.forte.pred <- rf.h.forte$test$predicted
     |      Out-of-bag   |       Test set    |
Tree |      MSE  %Var(y) |      MSE  %Var(y) |
  50 |    30.27    25.05 |    40.52    25.55 |
 100 |    29.12    24.10 |    39.43    24.87 |
 150 |    28.71    23.76 |    39.63    24.99 |
 200 |    28.46    23.56 |    39.86    25.14 |
 250 |    28.58    23.65 |    40.33    25.43 |
 300 |    29.03    24.03 |     40.8    25.73 |
 350 |    28.98    23.98 |    40.56    25.58 |
 400 |    28.71    23.76 |    40.65    25.63 |
 450 |    28.81    23.84 |    40.38    25.47 |
 500 |    28.48    23.57 |    40.63    25.62 |
In [108]:
class.forte.pred<- classify(rf.h.forte.pred)
y.forte.class <- classify(y_test_bin_Forte)
table.forte = table(pred.reg = class.forte.pred , observations = y.forte.class); table.forte
        observations
pred.reg  B  C  D  E  F  G
       D  3  2  8  0  0  0
       E  0  1 11 16  8  0
       F  0  0  0 12 19 11
       G  0  0  0  0  6 17

Modele de prediction : faibles consommations

In [109]:
rf.h.faible <- randomForest(x = x_train_bin_Faible, y = y_train_bin_Faible ,
                      xtest = x_test_bin_Faible, ytest = y_test_bin_Faible,
                     ntree = 500,do.trace = 50, importance = TRUE)

rf.h.faible.pred <- rf.h.faible$test$predicted
     |      Out-of-bag   |       Test set    |
Tree |      MSE  %Var(y) |      MSE  %Var(y) |
  50 |    17.82    44.73 |    19.22    51.95 |
 100 |    17.25    43.27 |    18.98    51.31 |
 150 |    17.06    42.82 |    18.51    50.04 |
 200 |    16.83    42.22 |    18.21    49.23 |
 250 |    16.76    42.05 |    18.01    48.68 |
 300 |    16.87    42.32 |    17.96    48.54 |
 350 |    16.84    42.25 |     18.1    48.93 |
 400 |    16.83    42.24 |     18.1    48.93 |
 450 |    16.82    42.20 |    18.12    48.98 |
 500 |    16.82    42.20 |    18.12    48.98 |
In [110]:
class.faible.pred<- classify(rf.h.faible.pred)
y.faible.class <- classify(y_test_bin_Faible)
table.faible = table(pred.reg = class.faible.pred , observations = y.faible.class); table.faible
        observations
pred.reg  A  B  C  D
       A 51 15  2  0
       B 10 19  6  0
       C  0  2 11  1

3. Mise en commun des modèles

In [111]:
pred.hierarc = c(class.faible.pred,class.forte.pred)
real.hierarc = c(y.faible.class,y.forte.class )
table.hierarc = table(pred.reg = pred.hierarc, observations = real.hierarc); table.hierarc
        observations
pred.reg  A  B  C  D  E  F  G
       A 51 15  2  0  0  0  0
       B 10 19  6  0  0  0  0
       C  0  2 11  1  0  0  0
       D  0  3  2  8  0  0  0
       E  0  0  1 11 16  8  0
       F  0  0  0  0 12 19 11
       G  0  0  0  0  0  6 17

4. Resultats

In [112]:
cat("Score : ",sum(diag(table.hierarc))/sum(table.hierarc), "\n" )
cat("Mean Absolute Error : ", MAE(table.hierarc), "\n")
cat("Macro F1 : ", mean(macro_F1(table.hierarc)) )
Score :  0.6103896 
Mean Absolute Error :  0.4155844 
Macro F1 :  0.5819292

Bilan : Les scores obtenus sont bons, juste en dessous de notre meilleur prédicteur des parties précédentes qui était le random Forest.

Premièrement, le modèle de régression logistique binaire pour prédire entre les consommations d'energie fortes et faibles est très bon (précision de 0.96).

En regardant la matrice de confusion, on voit que cette méthode de classification deux modèles qui apprennent spécifiquement sur des données de consommation faible ou fortes, n'apporte pas un tel gain car on se trompe encore beaucoup pour les classes intermédiaires (B,C, D et E). Plus particulièrment, on voit que la prédiction pour un individu se situant réellement dans la classe B est toujours compliquée, on va se tromper 2 fois sur 3 de prédiction en prédisant soit la classe A ou C.

On peut faire le lien avec ce que nous avons dit précedemment dans la conclusion, ces erreurs de prédictions sont dues aux intervalles choisies pour appartenir aux classes.

On peut finalement se demander si un tel modèle ne serait pas plus performant si on avait plus de données. Il se pourrait en effet que les modèles basés sur les consommations fortes et faibles soient un peu plus performants avec plus d'apprentissage.